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ABSTRACT 

Rotating magnetized compact objects and their accretion discs can generate strong 
toroidal magnetic fields driving highly magnetized plasmas into relativistic jets. Of significant 
concern, however, has been that a strong toroidal field in the jet should be highly unstable to 
the non-axisymmetric helical kink (screw) m - 1 mode leading to rapid disruption. In addi- 
tion, a recent concern has been that the jet formation process itself may be unstable due to 
the accretion of non-dipolar magnetic fields. We describe large-scale fully three-dimensional 
global general relativistic magnetohydrodynamic simulations of rapidly rotating, accreting 
black holes producing jets. We study both the stability of the jet as it propagates and the sta- 
biUty of the jet formation process during accretion of dipolar and quadrupolar fields. For our 
dipolar model, despite strong non-axisymmetric disc turbulence, the jet reaches Lorentz fac- 
tors of r ~ 10 with opening half-angle Oj ~ 5° at 10^ gravitational radii without significant 
disruption or dissipation with only mild substructure dominated by the m - \ mode. On the 
contrary, our quadrupolar model does not produce a steady relativistic (F > 3) jet due to mass- 
loading of the polar regions caused by unstable polar fields. Thus, if produced, relativistic jets 
are roughly stable structures and may reach up to an external shock with strong magnetic 
fields. We discuss the astrophysical implications of the accreted magnetic geometry playing 
such a significant role in relativistic jet formation, and we outline avenues for future work. 

Key words: accretion discs, black hole physics, galaxies: jets, gamma rays: bursts, MHD, 
instabilities, relativity, methods: numerical 



1 INTRODUCTION 

Astrophysical jets were discovered by Heber Curtis in 1917, 
, who described M87's jet as "a curious straight ray ... con- 
nected with the nucleus" jCurtij il918). The M87 jet is the 
\ most well-studie d of all Jets associated with active galactic nu- 
' clei (AGN) (e.g. iBaade & Minkowskilll954 iBlandford & Konigll 
ll979al : lR5ser & Meisenheimej|l993h . M87's jet structure has been 
observed down t o tens of gravitational r adii of the puta ti ve black 
hole (BH) (e.g. ljunoretal . 1999; Kov alev et al.l |2007|; ILv et al.l 
l2007h with impressive animations created i Walker et alj " 20081^ 
Jets have now also be en observed in many other AGN/blazars 
i Bridle&Pe"rievlll984). in neutron star and B H x-ray binaries 
I Mirabel & Rodriguez " 19991; iFender et alj|2004h . in Herbig-Haro 
objects (e.g. Konigl J982j) , and are required for gamma-ray bursts 
(GRBs) (e.g. lPiradl200^ . Challenges include explaining the jet 
formation process, the stability of jet formation and jet propaga- 
tion, how jets accelerate and collimate, and how jets obtain their 
composition and substructure both near and far from the central 
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object. Jet studies are complicated by the system's environment, 
such as how a jet must drill through a massive envelope in the col- 
lapsar model, while FRII jets extend up to hot spots. For quasar 
syst ems, jets play an impo rtant role in limiting BH mass growth 
(e.g.l Pi Matteo et alj2o6a) and d riving hot bubbles that limit cool- 
ing flows (e.g. McNamara et al.i r2005). However, the efficiency of 
the energy-momentum transfer remains unknown and probably de- 
pends on jet structure and stability. 

The most universally applicable jet paradigm involves some 
form of magnetic-driving with strong toroidal fields that form, 
accelerate, and (internally) collimate jets via magnetized accre- 
tion discs iLvnden-Belllll969l; IBlandford & Reeslll974 lLovelac3 
'1976'; 'Blandford & Pavnel Il982r or accreting, rotating BHs 
(Blandford & Znajek 1977, BZ). This paradigm was bolstered by 
the realization that accretion can be driven by magnetoro tational 
turbu lence generating a strong magnetic field tBalbus & HawlevI 
Il998l) . Especially if discs are thick near BHs, then stronger jets and 
winds are driven by either stronger turbulent magnetic fields for the 
other wise same mass accretion rate (e.g. Meier 2001; Mill er et al.l 
I2OO6I) or by large-scale fields adverted from large radii (see §3 & 
§4 in iLivio et alJI 199^ . A magnetic field may preserve jet com- 
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position against entrainment jRosen et" ai] |l999l) . The variations in 
magnetic field strength! and B H spin may explain the diversity of 
jet systems like FRI/FRIFs (JV[eiej|l999ll . although rotation mea- 
sures imply unexpected field orientations feavala & TavlorlllOOSi) 
and s imple models of decelerating jets fit FRIs (e.g. Laing et al] 



l2006i) . Soltan efficiency (and other) argu ments suggest quasars con- 
tain BHs that are rapidly spinnin g (e.g. Gammie et alj|20()4l) and 



maybe maximally spinning (e.g. (Allen et al.ll2006h . The intrinsic 
interest (and cosmological importance) of jets motivates testing 
whether the magnetic paradigm can explain their observed struc- 
ture and stability. 

Now roughly 90 years after Heber Curtis's discovery, the 
straightness of many observed jets remains as their most inex- 
plicable feature given, e.g., fusion devices show strong toroidal 
fields are viole ntly unstable to helical kink (screw) modes (e.g. 
lBatemanlll978l) . Astrophysic al jet stability rese arch has revealed 
a large number of modes (e.g. iKadomtsevll 19661) that can be unsta- 
ble including "reflection" resonant modes ('Pavne&Cohnl ll985b . 
Kelvin-Helmholtz (KH) modes (e.g. Ferrari et al. 1978 and refer- 
ences therein), and current-driven modes iBenfordll98ll) . With per- 
turbations of the form ^•(""t'+K+ii^-^t) ^ ^ nearly universal result from 
these simplified models is that the m = I kink mode is the most 
dangerous mode that could result in complete disruption and dissi- 
pation. 

Even if simplified jet models are kink mode unstable, they 
may b e stabilized by introducing gradual shear (e.g. i Mizuno et al.l 
20071 and references therein ), an external wind jHar dee & Hughes 
2003 *). sideways expansion (Rosen & Hardee 2000), and relativis- 
tic bulk motion. For some AGN jets, obser vations support a lack o f 
significant dissipation during propagation jSambruna et alj|2006h . 
If unstable, however, jets can be a source of heating, radiation, 
and hi gh-energy particles d ue to shocks (e.g. ^l andford & Konigl 



1979ah, recormection (e . g. Drenkhahn & Spruit! 120021 : ILvutikovi 



2006l :J Giannios & SDruilll2006h . viscous shear, turbulent cascade 
(e.g. iBegelmarT 1998t) . and a break-down of the ideal single- 



component fluid approximation l lTrussoni et al.|[l988h . 

For magnetized jets, the current-driven screw (n > 0,m = 1) 
mode is potentially most disruptive. For cylindrical force-free equi- 
libria one obtains the Kruskal-Shafranov (KS) instability criterion 
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where B* and are the toroidal and poloidal field strengths, 
R = r sin is cylindrical radius, and r is poloidal extent. This sug- 
gests jets are unstable beyond the Alfven surface where > B'' 
and r > R, located at only r < lOM (in this Letter, G = c = 1) for 
rotating BHs or accretion discs. The KS criterion implies jets go 
unstable before accel erating to relativis tic (F > 3) speeds as likely 
only after r ~ lOOM ( lMcKinnevl2006 al). and the KS criterion prob- 
ably cannot explain some FRIIs extending to r ~ lO^M. 

Advanced linear stability analyses from normal mode and 
extremal energy arguments for simplified cylindrical jets are 
often based upon restrictive assumptions, which ha s lead re- 



searchers to suggest th at jets a re violently unstable Begelmar 
19981 : iLvubarskiill 19991 : lLi H200(]|). mildly unstable ( e.g. 



200C), or even stable (e.g. Istomin & Parievlll99g) to the screw 
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mode. iTomimatsu et alj (|2001) found that relativi sti c field rotation 
of fre ely expanding solutions ( lMcKinnevll2006al l3: iNaravan et al.l 
l2007h allows jets to be unstable only if both the KS criterion and 
their criterion. 
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are satisfied, where Qf is the field line rotation frequency and c 
is the speed of light. This implies jets are marginally stable until 
a strong external medium interaction. Their analysis is suggestive, 
but it only strictly applies inside, not through, the Alfven surface. 
So far, no sufficiently general analytical screw stability analysis has 
been performed for magnetically-dominated relativistic jets. 

Analytical approaches become intractable for more realistic 
jets. It remains diffi cult to compare theory with obser vations (e.g. 
IWorrall et ^ l2007h and laboratory experiments (e.g. ICiardi et aU 
2008"). Primarily, numerical magnetohydrodynamical (MHD) sim- 
ulations have proven useful to study realistic jet models. Simu- 
latio ns range from injecting a n arbitrary jet from a surface inlet 
(e.g. 'Nakamura & Meier"2004': 'Zhang et al."2004j; iLeismann et"al] 
|2005: Komissarov et al. 2008; Moll et al. 20 0^ to injecting a 



jet fr om an unresolved Keplerian disc (e.g., ITchekhovskov et al 



2001 



2008 ) , and to evolving both the disc and jet (e.g., IMawtey et al. 
McKinney & Gammiell2002l. |2004 iKigure & Shibatall2005l: 

i.^^li. Itt„..,_. Tt.„ ,j k,_T...._..i yLr^A 
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McKinnev' '2005'; 'Hawlev & Krolik' '200^: iMcKinneV '2006d 
McKinnev & Naravan 2007a, b; Komissarov & McKinnev 2007]). 
Advanced 3D MHD simulations that inject jets from an inlet find 
that KH kink mode s are stabilized by sheaths around the jets 
jMizuno et al.l 120(371) and that even non-relativ istic screw mode s 
saturate before causing magnetic dissipation ("Moll et alj l2008h . 
More realistic simulations are crucial since analytical experience 
suggests free parameters in jet-injection simulations probably play 
a significant role. In particular, only global simulations allow a sta- 
bility study of the actual jet formation process in the presence of 
disc turbulence and different global field geometries. Accretion of 
small quadrupolar field loops was already shown to degrade the jet 
jMcKinnev & Gammiell2004l : lBeckwith et al.ll2008l) . but this could 
be due to their choice of starting with small field loops in the disc 
with numerical dissipation not allowing the development of a large- 
scale quadrupolar field. 



2 NUMERICAL MODEL 

We perform fully 3D global general relativistic MHD (GRMHD) 
simulations starting with an equilibrium matter torus, whose an- 
gular momentum is aligned with the BH (Kerr met ric) spin. To 
facilitate comparisons, we follow iMcKinnevI ( l2006a[) and choose 
a torus pressure maximum at r = 12M, inner edge at r = 6M, 
and adiabatic index y = 4/3 giving disc thickness 69 ~ ±0.3. 
For BH spins of a/M > 0.4, simulations o f such tori are qualita- 
tively similar jMcKirmev & Gammiell2004h . We choose all models 
to have a/M = 0.92 (hole angular frequency, Q.h = a/{2Mr+) » 
0.33M with horizon radius, r+) such that the BH i s in spin 
equilibrium for our disc thickness (Gammie et al. '2004"). We use 
the co nservative unsplit 3D GRMHD code HARM ( Ga mmie et al.l 
l2003l) . Kerr-Schild coordinates, 4th-o rder interpolation and 4th- 
order Runge-Kutta (lMcKinnevll2006al) . a robust inversion scheme 
dMignone & McKinnev 2007h . a staggered fi eld scheme (McK- 
inney et al., in prep.), a nd other advances dMcKinnevI [2006bl: 
ITchekhovskov et al.ll2007l) . 

We consider both dipolar and quadrupolar field geometries. 
The dipolar model starts with a single field loop within the torus as 
in McKinney ( 2006a). Dipolar models correspond to the most-often 
simulated jet (or jet+disc) model were the current sheet is assumed 
to be at (or develops near) the equator. To give the quadrupole ge- 
ometry the best chance of producing a jet, we study a large-scale 
quadrupolar field with vector potential (/> component 
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using a paraboloidal-like potential given by 

Ad,poie = (l/2)[(r + raff- + 2MU(i - ln(A))], (4) 

where /_ = 1 - cos'' 6, f+ = 1 + cos'' 6, v = 3/4, = 4, = 4, 
and applies for d < n/2 and for 9 > n/2 when letting d ^ n - d.ln 
this model, current sheets form above and below the equator. From 
prior GRMHD simulations, we expect primarily the initial field's 
multipole order to be important, and particular model parameter 
values should be unimportant once a quasi-steady state is reached. 
All models have initial gas pressure per unit magnetic pressure of 
w 100 at the equator in the disc. We allow the comoving magnetic 
energy per rest-mass energy up to only 100 during mass evacuation 
near the BH (see floor model in iMcKinnevll 2006ah . 

Spherical polar, not Cartesian, coordinates are used since pre- 
ferred for rotating jets. Our fiducial models have resolutio n 256 x 
128x32 in rxdX(p, with non-uniform grid as in McKinney l l2006ah . 
except Rq = and n,. = 1 in their equation (18). Based upon code 
tests, our 2nd-order monotonized central limiter scheme would re- 
quire roughly 4x the per-dimension resolution to obtain the accu- 
racy our 4th-order scheme by the end of the simulation. Unlike 
prior 3D GRMHD simulations, the grid warps to resolve the disc 
at small radii and follows the coUimating jet at large radii giving 
roughly 3x more angular resolution at large radii. Hence, compared 
to any scheme similar to the original 2nd-order HARM scheme, our 
effective resolution is roug hly 1024 x 1536 x 128 . Unlike most 3D 
GRMHD simulations (e.g. iBeckwith et alj|2008h . we include the 
full A(p = In extent as required to resolve the m = 1 mode and in- 
clude the full A6 = n extent (no cut-out at poles). As [Fragile et al.l 
we use transmissive (not reflecting) polar boundary condi- 
tions. As they state, the singularity need not be treated specially 
for centered quantities in a finite-volume scheme. Our field is stag- 
gered, and the polar value of 6^ is evolved by using the analytical 
limit of the finite volume induction equation at the pole such that 
angular-dependent area factors cancel (McKinney et al., in prep.). 
Coordinate directions twist at the pole leading to some dissipation, 
but this is significantly reduced by our 4th-order scheme that well- 
resolves up to m = 4 with 32 (p cells. At the iimer torus edge, cells 
have aspect ratio 1:5:10 and the fastest-gro wing magnetorotat ional 
mode is resolved with 6 cells, as sufficient jShafee et alj200^ . We 
also studied resolutions of 128 x 128 x 32, 128 x 64 x 32, and 
128 X 64 X 16; the jet's Fourier m = 1, 2, 3 power is converged to 
20%. Using 128 angular cells and staggered field scheme were re- 
quired for MHD jet invariants to be conserved to < 10%, w hich is 
evidence of an accurate solution ( Tchekhovsko v et alj2008l) . 

Most disc+jet simulations do not evolve to large enough radii 
to resolve a highly relativistic jet. For magnetically-dominated 
paraboloidal jets, the maximum Lorentz factor at large radii is 

r.o.3(-) , (5) 

jTchekhovskov et^l2008h . We choose an outer box radius of 
lO^M as required to reach F ~ 10. All simulations ran a duration of 
5000M, which is 192 orbits at the inner-most stable circular orbit 
(ISCO) (risco ~ 2.2M) and 50 orbits at the initial inner torus edge. 
The accretion rate of mass (M), energy, and angular momentum are 
roughly constant with radius out to r ~ lOM by / ~ 3000M, indi- 
cating the disc has reached a quasi-steady state. The slow/contact 
modes for the jet move with v/c > 0.2, so the jet is beyond the box 
by f = 5000M. We report many results at t ~ 4000M since this is 
before the jet partially reflects off the outer box. 




Figure 1. For dipolar model, shows inner ±100M cubical region with BH, 
accretion disc (pressure, yellow isosurface), outer disc and wind (log rest- 
mass density, low green, high orange, volume rendering), relativistic jet 
(Lorentz factor of F < 4, low blue, high red, volume rendering), and mag- 
netic field lines (green) threading BH. Despite non-axisymmetric turbu- 
lence, polar magnetically-dominated jets are launched by the BZ effect. 



3 RESULTS 

The fi ducial dipole model is overa l l similar to prior 2D simula- 
tions jMcKinnev & Gammid I2OO4I : iMcKinnevI l2006ij) . The BH- 
driven polar jet survives in a non-dissipated state to large radii. 
Each polar, magnetically-dominated jet at r+, 10, 10^, lO^M has 
constant electromagnetic luminosity of Lj a O.OlMc", with only 
a small secular drop as F in creases. This value is similar to higher 
resolution 2D simulations jMcKinnev & Gammid l2004h . The to- 
tal (disc+jet+wind) electromagnetic output peaks at r » lOM, 
but disc power is dissipated so does not survive at large radii 
jMcKinnev & Naravanl l2007iJ). Figure l[T]l shows the inner ±100M 
cubical region and Figure j2j shows out to z = lO^M by f = 4000M. 
The figures show the disc wind and relativistic jet generated by the 
rotating BH and magnetized, turbulent accretion disc. The jet is 
roughly stable out to z = lO'M reaching F ~ 5 - 10. Figure ^ 
shows the kinked polar jet structure of the poloidal current, RB'/', ca- 
pable of driving screw instabilities. We measure the Fourier power 
within the jet region defined by magnetic energy per rest-mass en- 
ergy, averaged for all (p, greater than one. At large distances the 
ni = 1,2,3,4 powers relative to ra = are 7%, 1%, 0.7%, 0.6%c 
in magnetic energy, 6%, 4%, 0.5%, 0.2% in Lorentz factor, roughly 
37%),7%,3%>,4%) in both rest-mass density (po) and RB'^, and 
20%), 13%), 7%, 6% in internal energy density. Both po and RB^ 
reach m = 1 power of 100% in the jet next to the outer disc edge 
at r = 20M. There is no indication of growth beyond perturbations 
induced by the disc turbulence, which appears to be the primary 
origin of jet substructure. 

Now we discuss our fiducial large-scale quadrupole model. 
GRMHD simulations show that no strong jet emerges due to the 
accre tion of higher multipole moments put initially within t he 
disc jMcKinnev & Gammidl2004lMcKinnev & Naravai]|2007j|bl: 
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Figure 2. For dipolar model, shows accreting BH generating relativistic jet 
(only one side shown) vertically out to IO'M within 8 « ±20° (350M X 
350M) aX t = 4000A/. Shown are outskirts of disc and wind (log internal 
energy density, cyan volume rendering), outer/inner boundary of perturbed 
jet and fragments of disc wind {RB'^, blue isosurface), and relativistic jet 
(Lorentz factor of F g 10, orange volume rendering) collimated within half- 
angle 0j a 5°. Despite perturbations, the jet is a stable structure. 



iBeckwith et"aDl2008h . In our fully 3D simulations, even putting in 
an initial large-scale quadrupolar field leads to no strong jet once 
equatorial symmetry is broken by / ~ 2500M. While the polar field 
strength is similar to that in the dipole model when the field threads 
the horizon, the polar regions are mass-loaded when tearing coronal 
current sheets eject polar field to slightly larger radii. Then, mag- 
netic pressure no longer balances against the low angular momen- 
tum disc material that moves into the polar region. The coronal- 
polar field is accreted and ejected throughout the simulation, which 
leaves no time allowed for the funnel to drain. This leads to order 
unity magnetic energy and internal energy per particle rest-mass 
energy, which is insuf ficient to generat e a high ly relativistic (F > 3) 
jet. Similar results by IBeckwith et alj 1200S) using a non-energy- 
conserving code suggest detailed thermodynamics do not control 
this process. At late times, (total and polar) BH electromagnetic 
power are negative. The polar regions have both inflows and out- 
flows, and there is only a disc-mass-loaded wind with an electro- 
magnetic power output per pole at large radii of L„. a 0.002Mc-, 
which is significantly less powerful than the dipole model. Similar 
as the disc wind in the dipolar model, the outflow has a weak disor- 
ganized poloidal field and a more organized toroidal field stronger 
by factors of typically 10 - 40 both near the BH and at large radii. 
Corresponding 2D simulations show less drastic, but comparable, 
degradation of the jet. A thinner disc may not allow as much mass- 
transfer to the poles, but thinner discs have weaker turbulent fields 
and inward advection of strong ordered field may not be possible 



for thin discs. Also, higher resolutions may lead to less vigorous 
reconnection or may show a more narrow, polar jet still emerges. 



4 DISCUSSION 

We have performed fully 3D global GRMHD simulations of accret- 
ing, rapidly rotating BHs and found that dipolar fields near BHs 
can launch magnetically-dominated, relativistic (F > 3) jets that 
survive to lO^M without significant disruption or measurable dis- 
sipation. Disc turbulence appears to be the primary cause of jet 
substructure that is dominated by the m = 1 mode, which has no 
measurable growth within the jet. Prior work applying a form of the 
Kruskal-Shafranov criterion (solution for non-relativistic , cylindri- 
cal equilibria) to highly magne tized relativistic flows (e.g. lLvutikovl 
|2006| -[g larmios & SDruitl200^ . needs to be reevaluated to consider 
the stabilizing effects of field rotation, gradual shear, a surrounding 
sheath, and sideways expansion as present in the simulations. Un- 
like dipolar fields, quadrupolar fields near BHs lead to only weak, 
turbulent outflows and negligible magnetically-dominated polar re- 
gions and no relativistic (F > 3) jets. Since our simulations with 
relativistic jets have no current sheets within the jet, reconnection 
may not be an important source of dissipation unlike assumed by 
some models (e.g. Drenkhahn & Spruit.2002 .) 

These and prior GRMHD simulation results suggest that a ro- 
tating (a/M > 0.4) BH is a necessary, but not sufficient, condition 
to produce a highly relativistic (F > 3) jet. In addition, one re- 
quires the accreted magnetic field to be mostly dipolar, rather than 
hig her-order, so a dipola r field threads the region near the BH (see 
also lNaravan et al .120031) . This might explain various observations, 
such as the dichotomy of FRI and FRII systems. FRI's are found 
in rich clusters, are two-sided so weakly relativistic, and have dis- 
sipative emission near the core. FRII's are found in poor groups 
or isolated, are one-sided so more re lativistic, are more pow erful, 
and dissipate little till the radio lobe ( IOwen&Led"iowlll994h . The 
FRI/FRII dichotomy may then be due to the complexity of the en- 
vironment (e.g. through hierarchical merging) controlling the field 
multipole structure. Then, FRII systems are primarily BH-driven 
able to pierce through an ambient medium, while FRI systems 
are those mostly driven by the broader, dissipative, magnetically- 
disordered disc wind with F < 3 that one expects to be more easily 
entrained, slowed , and disrupted, as consistent with observations 
jLaing et al.ll2006l) . Radial structure (e.g. arcs and knots) could be 
due to accretion switching between dipolar and higher-order multi- 
poles. For M87, there could be a dark or boosted relativistic spine 
with the slower, dissip ative disc wind p roducing emission on scales 
within several parsecs jKovalev et alBoO?,) . For SrgA*, no jet may 
emerge because of accretion from various stellar clusters generat- 
ing a dominant non-dipolar field jNavakshin et al.ir2007h . For x- 
ray binaries, jets in the low-hard states could be driven by dipolar 
fields that could even accumulate to the point of lowering accretion 
rates jlgumenshchev et aill2003h , intermediate to soft states could 
involve higher-order multipole moments, and transient jets from the 
hard-to-soft transitions could occur due to dissipation of the dipolar 
component. For GRBs, the BH-disc system may be required to be 
highly symmetric to maintain a strong dipolar field to produce an 
ultrarelativistic jet. That ordered poloidal field must be accreted as- 
sumes no dynamo exists for generati ng a baryon-pure, larg e-scale 
poloidal field from disorganized field IBeckwith et al.ll2008l) . 

Future jet studies should consider the eflFects of much higher 
resolutions, misaligned BH-disc accretion (since misaligned sys- 
tems may more readily produce non-dipolar fields), larger radii of 
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lO'M for AGN and lO'-M for GRBs (to determine very large-scale 
stability and to obtain larger T), resistivity and viscosity, disc radial 
extent (that limits the terminal Lorentz factor since the lack of the 
disc and supportive disc wind allows the jet to become monopo- 
lar and so accelerate inefficiently), disc thickness (that can control 
the strength of turbulent or advected field), other magnetic field 
geometries (including with non-zero net helicity), BH spin (espe- 
cially very low and very high), cooling (such as neutrino cooling in 
collapsar discs), and the presence of an extended massive envelope 
as in the collapsar model (freely expanding outflows simulated here 
apply to a late phase after the jet drills through the envelope). Future 
studies should also do a quantitative analysis of the modes within 
the jet to identify which mode types are present. The simulated jet 
can be used as a well-motivated background state for future lin- 
ear perturbation analyses, parameter searches, and synchrotron and 
inverse Compton maps for, e.g., VLBI, Chandra, and Fermi. 
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